##########################################################################
## Replication syntax for Hegewald (2023, JEPP):
## "Locality as a safe haven: place-based resentment and political trust in local and national institutions"
##########################################################################

# R version: 4.2.2 
# Platform: x86_64-w64-mingw32/x64 (64-bit)

# 1. Load packages #######################################################
library(here)
library(foreign)
library(tidyverse)
library(magrittr)
library(psych)
library(reshape2)
library(ggpubr)
library(modelsummary)
library(interactions)
library(stargazer)

# Package version info: 
# stargazer_5.2.3    interactions_1.1.5 modelsummary_1.4.2 ggpubr_0.5.0       reshape2_1.4.4     psych_2.2.9        foreign_0.8-83    
# magrittr_2.0.3     forcats_0.5.2      stringr_1.5.0      dplyr_1.0.10       purrr_0.3.5        readr_2.1.3        tidyr_1.2.1       
# tibble_3.1.8       ggplot2_3.4.0      tidyverse_1.3.2    here_1.0.1  

# 2. Basic setup #########################################################
rm(list=ls())
i_am("Hegewald_JEPP_replication.R")

# 3. Import data ##########################################################
dat_full <- read.csv((here("dat_full.csv")), sep=",") ### load in survey data
LAI <- read.spss((here("Tot_File_1990_2020_Temp5.sav")), to.data.frame = T) ### load in LAI data from Ladner et al. (2021) 

# 4. Prepare data & variables #############################################
# 4.1. Drop non-binary respondents and create respondent id ###############
dat_full %<>% filter(gender!="non-binary") ### 11 observations dropped
dat_full$pers_id <- seq.int(nrow(dat_full)) ### create respondent id

# 4.2. Create local trust differential ####################################
dat_full$trstdiff <- dat_full$loctrst-dat_full$nattrst 

# 4.3. Create Place-based resentment scale ################################
psych::alpha(dat_full[,c(13:17)]) ### 0.83 

pbr.fa <- fa(dat_full[,c(13:17)], nfactors = 1)

colnames(pbr.fa$loadings) <- c("Place-based resentment")
rownames(pbr.fa$loadings) <- c("Economic", "Representation (A)", "Representation (B)", "Culture (A)", "Culture (B)")

fa.diagram(pbr.fa, digits = 2, main = "")

dat_full$resent_munis.std <- as.numeric(pbr.fa$scores)

summary(dat_full$resent_munis.std)
sd(dat_full$resent_munis.std)

dat_full$resent_munis.std <- as.numeric(scale(dat_full$resent_munis.std)) 

summary(dat_full$resent_munis.std)
sd(dat_full$resent_munis.std)

# 4.4. Create component indicators ########################################
dat_full$eco.std <- as.numeric(scale(dat_full$eco))

psych::alpha(dat_full[,c(14:15)]) ### 0.81 

dat_full$rep_av <- (dat_full$rep1 + dat_full$rep2)/2
dat_full$rep.std <- as.numeric(scale(dat_full$rep_av)) 

summary(dat_full$rep.std)
sd(dat_full$rep.std)

psych::alpha(dat_full[,c(16:17)]) ### 0.67 

dat_full$cult_av <- (dat_full$cult1 + dat_full$cult2)/2
dat_full$cult.std <- as.numeric(scale(dat_full$cult_av)) 

summary(dat_full$cult.std)
sd(dat_full$cult.std)

# 4.4. Standardize control variables ######################################
dat_full$age.std <- as.numeric(scale(dat_full$age))
dat_full$lr.std <- as.numeric(scale(dat_full$lr))
dat_full$soctrst.std <- as.numeric(scale(dat_full$soctrst))
dat_full$plid.std <- as.numeric(scale(dat_full$plid))

# 4.5. Set base of urban-rural indicators #################################
dat_full$rur <- ifelse(dat_full$rur=="urban", "Urban residence", "Rural residence")
dat_full$rur <- as.factor(dat_full$rur)
dat_full <- within(dat_full, rur <- relevel(rur, ref = "Urban residence"))

dat_full$rur_cont <- as.factor(dat_full$rur_cont)
dat_full$rur_cont <- ordered(dat_full$rur_cont, levels = c("very urban", "rather urban", "rather rural", "very rural"))
levels(dat_full$rur_cont) <- c("Very urban", "Rather urban", "Rather rural", "Very rural")

dat_full$rur_cont <-  factor(dat_full$rur_cont, ordered = FALSE)
dat_full <- within(dat_full, rur_cont <- relevel(rur_cont, ref = "Very urban"))

dat_full$peri <- as.factor(dat_full$peri)
dat_full <- within(dat_full, peri <- relevel(peri, ref = "center"))

# 4.6. Set base of gender, educ and opposition ############################
dat_full$gender <- as.factor(dat_full$gender)
dat_full <- within(dat_full, gender <- relevel(gender, ref = "male"))

dat_full$educ <- as.factor(dat_full$educ)
dat_full <- within(dat_full, educ <- relevel(educ, ref = "Low"))

dat_full$vote_opp <- as.factor(dat_full$vote_opp)
dat_full <- within(dat_full, vote_opp <- relevel(vote_opp, ref = "government"))

# 5. Define themes ########################################################
theme_own1 <- theme_minimal() +
  theme(text=element_text(size=20, colour="black"),
        axis.title.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.x = element_text(color="black", 
                                   size=20),
        axis.text.y = element_text(colour="black"),
        legend.pos="none")

theme_own2 <- theme_minimal() +
  theme(text=element_text(size=25, colour="black"),
        axis.text = element_text(colour="black"),
        axis.ticks = element_line(colour="black"),
        legend.title = element_blank(),
        legend.pos="bottom")

# 6. Main Analysis  #######################################################
# 6.1 Figure 1 ############################################################
t.test(dat_full$loctrst, dat_full$nattrst, paired = TRUE, alternative = "two.sided") 

dat_trst <- dat_full %>% select(pers_id, nattrst, loctrst)
dat_trst_long <- melt(dat_trst, id=c("pers_id")) 
dat_trst_long$variable <- ifelse(dat_trst_long$variable=="nattrst", "National trust", "Local trust")

Fig1A <- ggplot(dat_trst_long, aes(variable, value)) +
  stat_boxplot(aes(variable, value), width=0.5)+  
  geom_boxplot(aes(variable, value),outlier.shape=1, width=0.5) +    
  stat_summary(fun=mean, geom="point", size=2, color="#FF0000") + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", color = "#FF0000", width = 0.1) + 
  scale_color_manual(values = c("#666666", "#000000")) +
  stat_compare_means(label = "p.format", 
                     method = "t.test", paired = TRUE, comparisons = list(c("National trust", "Local trust"))) + 
  labs(y = "Local trust / National trust", title = "A: Full sample") + theme_own1
Fig1A

dat_trst_long %>%
  group_by(variable) %>%
  summarise(mean = mean(value))

Fig1B_pre <- dat_full %>% 
  ggplot(aes(trstdiff)) + 
  geom_histogram()  + 
  geom_histogram(binwidth = 1)+labs(title="B: Full sample",x="Local trust differential", y = "Frequency") + theme_own2

trst.diff_mean <- dat_full %>% summarise(grp.mean=mean(trstdiff))

Fig1B <- Fig1B_pre+geom_vline(data=trst.diff_mean, aes(xintercept=grp.mean),
                              linetype="dashed", color = "red", linewidth=1)
Fig1B

dat_trst <- dat_full %>% select(pers_id, rur, nattrst, loctrst)
dat_trst %<>% filter(rur=="Urban residence")
dat_trst$rur <- NULL

summary(dat_trst$nattrst)
summary(dat_trst$loctrst)

t.test(dat_trst$loctrst, dat_trst$nattrst, paired = TRUE, alternative = "two.sided") 

dat_trst_long <- melt(dat_trst, id=c("pers_id")) 
dat_trst_long$variable <- ifelse(dat_trst_long$variable=="nattrst", "National trust", "Local trust")

Fig1C <- ggplot(dat_trst_long, aes(variable, value)) +
  stat_boxplot( aes(variable, value), width=0.5)+  
  geom_boxplot( aes(variable, value),outlier.shape=1, width=0.5) +    
  stat_summary(fun=mean, geom="point", size=2, color="#FF0000") + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", color = "#FF0000", width = 0.1) + 
  scale_color_manual(values = c("#666666", "#000000")) +
  stat_compare_means(label = "p.format", 
                     method = "t.test", paired = TRUE, comparisons = list(c("National trust", "Local trust"))) + 
  labs(y = "Local trust / National trust", title = "C: Urban residence") + theme_own1
Fig1C

dat_trst <- dat_full %>% select(pers_id, rur, nattrst, loctrst)
dat_trst %<>% filter(rur=="Rural residence")
dat_trst$rur <- NULL

summary(dat_trst$nattrst)
summary(dat_trst$loctrst)

t.test(dat_trst$loctrst, dat_trst$nattrst, paired = TRUE, alternative = "two.sided") 

dat_trst_long <- melt(dat_trst, id=c("pers_id")) 
dat_trst_long$variable <- ifelse(dat_trst_long$variable=="nattrst", "National trust", "Local trust")

Fig1D <- ggplot(dat_trst_long, aes(variable, value)) +
  stat_boxplot( aes(variable, value), width=0.5)+  
  geom_boxplot( aes(variable, value),outlier.shape=1, width=0.5) +    
  stat_summary(fun=mean, geom="point", size=2, color="#FF0000") + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", color = "#FF0000", width = 0.1) + 
  scale_color_manual(values = c("#666666", "#000000")) +
  stat_compare_means(label = "p.format", 
                     method = "t.test", paired = TRUE, comparisons = list(c("National trust", "Local trust"))) + 
  labs(y = "Local trust / National trust", title = "D: Rural residence") + theme_own1
Fig1D

ggarrange(Fig1A, Fig1B, Fig1C, Fig1D,
          ncol = 2, nrow = 2)

# 6.2 Table 3 #############################################################
FE_add = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes"), ncol = 4)) 
attr(FE_add, "position") = 21

Table3 <- list(
  "(1)" = lm(trstdiff ~ resent_munis.std + cntry, data = dat_full),
  "(2)" = lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + cntry, data = dat_full),
  "(3)" = lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = dat_full)
)

modelsummary(Table3, fmt = 3, stars = T,
             title = 'OLS regression results: Local trust differential on place-based resentment scale.', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "lr.std" = "Left-right (Std.)",
                          "soctrst.std" = "Social trust (Std.)",
                          "plid.std" = "Place-based identity (Std.)",
                          "vote_oppopposition" = "Opposition vote (b.=government vote)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.3 Figure 2 & Table A5 #################################################
int <- lm(trstdiff ~ resent_munis.std*rur + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = dat_full)

Fig2 <- interact_plot(model = int, pred = resent_munis.std, modx = rur, interval = TRUE, int.type = "confidence",
              int.width = 0.95, colors = "Greys") 

Fig2 <- Fig2 + geom_rug(data=dat_full, aes(x = resent_munis.std), inherit.aes = F) +
               labs(y = "Local trust differential",
               x = "Place-based resentment (Std.)") + theme_own2
Fig2

TableA5 <- list(
  "(1)" = lm(trstdiff ~ resent_munis.std*rur + cntry, data = dat_full),
  "(2)" = lm(trstdiff ~ resent_munis.std*rur + gender + age.std + educ + inc_dec + cntry, data = dat_full),
  "(3)" = lm(trstdiff ~ resent_munis.std*rur + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = dat_full)
)

attr(FE_add, "position") = 25
modelsummary(TableA5, fmt = 3, stars = T,
             title = 'OLS regression results: Local trust differential on place-based resentment scale conditional on urban-rural residence (binary variable).', 
             add_rows = FE_add,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "lr.std" = "Left-right (Std.)",
                          "soctrst.std" = "Social trust (Std.)",
                          "plid.std" = "Place-based identity (Std.)",
                          "vote_oppopposition" = "Opposition vote (b.=government vote)",
                          "resent_munis.std:rurRural residence" = "Place-based resentment (Std.) X Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.4 Figure 3 & Table A7 #################################################
comp_rur <- lm(trstdiff ~ eco.std + rep.std + cult.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
dat_comp_rur <- broom::tidy(comp_rur)
dat_comp_rur %<>% filter(term=="eco.std" | term=="rep.std" | term=="cult.std")
dat_comp_rur$sample <- "Rural residence"

comp_urb <- lm(trstdiff ~ eco.std + rep.std + cult.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = subset(dat_full, dat_full$rur=="Urban residence"))
dat_comp_urb <- broom::tidy(comp_urb)
dat_comp_urb %<>% filter(term=="eco.std" | term=="rep.std" | term=="cult.std")
dat_comp_urb$sample <- "Urban residence"

dat_comp_fin <- rbind(dat_comp_rur, dat_comp_urb)
dat_comp_fin$term[dat_comp_fin$term=="eco.std"] <- "Economic (Std.)"
dat_comp_fin$term[dat_comp_fin$term=="rep.std"] <- "Representation (Std.)"
dat_comp_fin$term[dat_comp_fin$term=="cult.std"] <- "Cultural (Std.)"

dat_comp_fin$sample <- as.factor(dat_comp_fin$sample)
dat_comp_fin$sample <- fct_rev(dat_comp_fin$sample)
dat_comp_fin$estimate <- as.numeric(dat_comp_fin$estimate)
dat_comp_fin$std.error <- as.numeric(dat_comp_fin$std.error)
dat_comp_fin$var <- dat_comp_fin$term
dat_comp_fin$var <- factor(as.factor(dat_comp_fin$var), levels = c("Cultural (Std.)", "Representation (Std.)", "Economic (Std.)"))

interval1 <- -qnorm((1-0.95)/2)  # 95% multiplier
interval2 <- -qnorm((1-0.99)/2)  # 99% multiplier

Fig3 <- ggplot(dat_comp_fin, aes(x = var, y = estimate))
Fig3 <- Fig3 + geom_hline(yintercept = 0, color = "red", lty = 2)

Fig3 <- Fig3 + geom_linerange(aes(x = var, ymin = estimate - std.error*interval1,
                                ymax = estimate + std.error*interval1),
                              linewidth = 1, position = position_dodge(width = 1/2))

Fig3 <- Fig3 + geom_pointrange(aes(x = var, y = estimate, ymin = estimate - std.error*interval2,
                                 ymax = estimate + std.error*interval2),
                               linewidth = 1/2, position = position_dodge(width = 1/2))


Fig3 <- Fig3 + xlab("") + ylab("Estimated coefficient") + coord_flip() + theme_bw(base_size = 20)
Fig3 <- Fig3 + facet_wrap(vars(sample))

Fig3

FE_add2 = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 7)) 
attr(FE_add2, "position") = 25

TableA7 <- list(
  "(1)" = lm(trstdiff ~ eco.std + rep.std + cult.std + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
  "(2)" = lm(trstdiff ~ eco.std + rep.std + cult.std + gender + age.std + educ + inc_dec + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
  "(3)" = lm(trstdiff ~ eco.std + rep.std + cult.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = subset(dat_full, dat_full$rur=="Rural residence")),
  "(4)" = lm(trstdiff ~ eco.std + rep.std + cult.std + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
  "(5)" = lm(trstdiff ~ eco.std + rep.std + cult.std + gender + age.std + educ + inc_dec + cntry, data = subset(dat_full, dat_full$rur=="Urban residence")),
  "(6)" = lm(trstdiff ~ eco.std + rep.std + cult.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = subset(dat_full, dat_full$rur=="Urban residence"))
)

modelsummary(TableA7, fmt = 3, stars = T,
             title = 'OLS regression results: Local trust differential on place-based resentment dimensions (rural sample versus urban sample).', 
             add_rows = FE_add2,
             coef_omit = "cntry", 
             coef_map = c("eco.std" = "Economic (Std.)",
                          "rep.std" = "Representation (Std.)",
                          "cult.std" = "Cultural (Std.)", 
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "lr.std" = "Left-right (Std.)",
                          "soctrst.std" = "Social trust (Std.)",
                          "plid.std" = "Place-based identity (Std.)",
                          "vote_oppopposition" = "Opposition vote (b.= government vote)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 6.5 Figure 4 & Table A8 #################################################
LowLAI <- dat_full %>% filter(cntry=="HU" | cntry=="CZ" | cntry=="PL" | cntry=="EL" | cntry=="IT")
HighLAI <- dat_full %>% filter(cntry=="DE" | cntry=="ES" | cntry=="DK" | cntry=="FR")

Low_urb <- lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = subset(LowLAI, LowLAI$rur=="Urban residence"))
dat_Low_urb <- broom::tidy(Low_urb)
dat_Low_urb %<>% filter(term=="resent_munis.std")
dat_Low_urb$sample1 <- "Lower local autonomy"
dat_Low_urb$sample2 <- "Urban residence"

Low_rur <- lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = subset(LowLAI, LowLAI$rur=="Rural residence"))
dat_Low_rur <- broom::tidy(Low_rur)
dat_Low_rur %<>% filter(term=="resent_munis.std")
dat_Low_rur$sample1 <- "Lower local autonomy"
dat_Low_rur$sample2 <- "Rural residence"

High_urb <- lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, subset(HighLAI, HighLAI$rur=="Urban residence"))
dat_High_urb_rur <- broom::tidy(High_urb)
dat_High_urb_rur %<>% filter(term=="resent_munis.std")
dat_High_urb_rur$sample1 <- "Higher local autonomy"
dat_High_urb_rur$sample2 <- "Urban residence"

High_rur <- lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, subset(HighLAI, HighLAI$rur=="Rural residence"))
dat_High_rur <- broom::tidy(High_rur)
dat_High_rur %<>% filter(term=="resent_munis.std")
dat_High_rur$sample1 <- "Higher local autonomy"
dat_High_rur$sample2 <- "Rural residence"

dat_LAI_fin1 <- rbind(dat_Low_urb, dat_Low_rur)
dat_LAI_fin2 <- rbind(dat_LAI_fin1, dat_High_urb_rur)
dat_LAI_fin3 <- rbind(dat_LAI_fin2, dat_High_rur)
dat_LAI_fin <- dat_LAI_fin3

dat_LAI_fin$term[dat_LAI_fin$term=="resent_munis.std"] <- "Place-based resentment (Std.)"

dat_LAI_fin$term <- as.factor(dat_LAI_fin$term)
dat_LAI_fin$sample1 <- as.factor(dat_LAI_fin$sample1)
dat_LAI_fin$sample1 <- factor(dat_LAI_fin$sample1, levels = c("Lower local autonomy", "Higher local autonomy"))
dat_LAI_fin$sample2 <- as.factor(dat_LAI_fin$sample2)
dat_LAI_fin$sample2 <- fct_rev(dat_LAI_fin$sample2)
dat_LAI_fin$estimate <- as.numeric(dat_LAI_fin$estimate)
dat_LAI_fin$std.error <- as.numeric(dat_LAI_fin$std.error)

Fig4 <- ggplot(dat_LAI_fin, aes(x = term, y = estimate))
Fig4 <- Fig4 + geom_hline(yintercept = 0, color = "red", lty = 2)

Fig4 <- Fig4 + geom_linerange(aes(x = term, ymin = estimate - std.error*interval1,
                                ymax = estimate + std.error*interval1, color = sample2),
                              linewidth = 1, position = position_dodge(width = -1/2)) +
  scale_color_manual(values = c("black", "red"))

Fig4 <- Fig4 + geom_pointrange(aes(x = term, y = estimate, ymin = estimate - std.error*interval2,
                                 ymax = estimate + std.error*interval2, color = sample2, shape = sample2),
                               linewidth = 1/2, position = position_dodge(width = -1/2)) +
  scale_color_manual(values = c("black", "red"))

Fig4 <- Fig4 + xlab("") + ylab("Estimated coefficient") + coord_flip() + theme_bw(base_size = 20)

Fig4 <- Fig4 + facet_wrap(vars(sample1)) + 
  theme(legend.position="bottom") +
  theme(legend.title=element_blank())

Fig4

TableA8 <- list(
  "(1)" = lm(trstdiff ~ resent_munis.std + cntry, data = subset(LowLAI, LowLAI$rur=="Urban residence")),
  "(2)" = lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + cntry, data = subset(LowLAI, LowLAI$rur=="Urban residence")),
  "(3)" = lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = subset(LowLAI, LowLAI$rur=="Urban residence")),
  "(4)" = lm(trstdiff ~ resent_munis.std + cntry, data = subset(LowLAI, LowLAI$rur=="Rural residence")),
  "(5)" = lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + cntry, data = subset(LowLAI, LowLAI$rur=="Rural residence")),
  "(6)" = lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = subset(LowLAI, LowLAI$rur=="Rural residence")),
  "(7)" = lm(trstdiff ~ resent_munis.std + cntry, subset(HighLAI, HighLAI$rur=="Urban residence")),
  "(8)" = lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + cntry, subset(HighLAI, HighLAI$rur=="Urban residence")),
  "(9)" = lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, subset(HighLAI, HighLAI$rur=="Urban residence")),
  "(10)" = lm(trstdiff ~ resent_munis.std + cntry, subset(HighLAI, HighLAI$rur=="Rural residence")),
  "(11)" = lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + cntry, subset(HighLAI, HighLAI$rur=="Rural residence")),
  "(12)" = lm(trstdiff ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, subset(HighLAI, HighLAI$rur=="Rural residence"))
)

FE_add3 = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 13)) 
attr(FE_add3, "position") = 21

modelsummary(TableA8, fmt = 3, stars = T,
             title = 'OLS regression results: Local trust differential on place-based resentment scale (urban vs rural; high versus low local autonomy).', 
             add_rows = FE_add3,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "lr.std" = "Left-right (Std.)",
                          "soctrst.std" = "Social trust (Std.)",
                          "plid.std" = "Place-based identity (Std.)",
                          "vote_oppopposition" = "Opposition vote (b.= government vote)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

# 7. Appendix #############################################################
# 7.1 Table A1 ############################################################
dat_summary <- dat_full %>% select(gender, age, age.std, educ, inc_dec, loctrst, nattrst, trstdiff,
                                   resent_munis.std, eco, eco.std, rep1, rep2, rep_av, rep.std, 
                                   cult1, cult2, cult_av, cult.std, lr, lr.std, soctrst, soctrst.std,
                                   plid, plid.std, vote_opp, rur, rur_cont, peri)

dat_summary$gender.num <- as.numeric(dat_summary$gender)
dat_summary$gender.num <- ifelse(dat_summary$gender.num==1, 0, 1)

dat_summary$educ.num <- as.numeric(dat_summary$educ)
dat_summary$educ.num <- ifelse(dat_summary$educ.num==1, 0, 1)

dat_summary$rur.num <- as.numeric(dat_summary$rur)
dat_summary$rur.num <- ifelse(dat_summary$rur.num==1, 0, 1)

dat_summary$rur_cont <- as.character(dat_summary$rur_cont)
dat_summary$rur_rec <- recode(dat_summary$rur_cont, "Very urban"="0", "Rather urban"="1",
                              "Rather rural"="2", "Very rural"="3")
dat_summary$rur_cont <- as.numeric(dat_summary$rur_rec)

dat_summary$vote_opp.num <- as.numeric(dat_summary$vote_opp)
dat_summary$vote_opp.num <- ifelse(dat_summary$vote_opp.num==1, 0, 1)

dat_summary$peri.num <- as.numeric(dat_summary$peri)
dat_summary$peri.num <- ifelse(dat_summary$peri.num==1, 0, 1)

stargazer(dat_summary, digits=2, type = "text")

# 7.2 Figure A2 & Figure A3 ###############################################
valid <- dat_full %>% select(eco, rep1, rep2, cult1, cult2, loctrst, nattrst)

validation <- fa(valid[,c(1:7)], nfactors = 1)
validation ### BIC = 3684.34
colnames(validation$loadings) <- c("Place-based resentment/\nPolitical trust")
rownames(validation$loadings) <- c("Economic", "Representation (A)", "Representation (B)", "Culture (A)", "Culture (B)",
                                   "Local trust", "National trust")

fa.diagram(validation, digits = 2, main = "", cut = 0)

validation <- fa(valid[,c(1:7)], nfactors = 2)
validation ### BIC = 394.73
colnames(validation$loadings) <- c("Place-based resentment", "Political trust")
rownames(validation$loadings) <- c("Economic", "Representation (A)", "Representation (B)", "Culture (A)", "Culture (B)",
                                   "Local trust", "National trust")

fa.diagram(validation, digits = 2, main = "")

# 7.3 Table A4 & Figure A4 #################################################
TableA4 <- list(
  "(1)" = lm(trstdiff ~ resent_munis.std*rur_cont + cntry, data = dat_full),
  "(2)" = lm(trstdiff ~ resent_munis.std*rur_cont + gender + age.std + educ + inc_dec + cntry, data = dat_full),
  "(3)" = lm(trstdiff ~ resent_munis.std*rur_cont + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = dat_full)
)

attr(FE_add, "position") = 33
modelsummary(TableA4, fmt = 3, stars = T,
             title = 'OLS regression results: Local trust differential on place-based resentment scale conditional on urban-rural residence (full variable).', 
             add_rows = FE_add,
             coef_omit = "cntry",
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "rur_contRather urban" = "Rather urban",
                          "rur_contRather rural" = "Rather rural",
                          "rur_contVery rural" = "Very rural",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "lr.std" = "Left-right (Std.)",
                          "soctrst.std" = "Social trust (Std.)",
                          "plid.std" = "Place-based identity (Std.)",
                          "vote_oppopposition" = "Opposition vote (b.= government vote)",
                          "resent_munis.std:rur_contRather urban" = "Place-based resentment (Std.) X Rather urban",
                          "resent_munis.std:rur_contRather rural" = "Place-based resentment (Std.) X Rather rural",
                          "resent_munis.std:rur_contVery rural" = "Place-based resentment (Std.) X Very rural",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))


int_fullcont <- lm(trstdiff ~ resent_munis.std*rur_cont + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = dat_full)

FigA4 <- interact_plot(model = int_fullcont, pred = resent_munis.std, modx = rur_cont, interval = TRUE, int.type = "confidence",
                      int.width = 0.95, colors = "Greys") 

FigA4 <- FigA4 + geom_rug(data=dat_full, aes(x = resent_munis.std), inherit.aes = F) +
  labs(y = "Local trust differential",
       x = "Place-based resentment (Std.)") + theme_own2

FigA4

# 7.4 Figure A5 ###########################################################
LAI %<>% select(country_code, Country_name_57, LAI_Index_D7w_2020s)
LAI_EU <- LAI %>% filter(country_code == 2 | country_code == 3 | country_code == 4 | country_code == 5 | country_code == 6 | 
                country_code == 7 | country_code == 8 | country_code == 9 | country_code == 10 | country_code == 11 | 
                country_code == 13 | country_code == 14 | country_code == 15 | country_code == 17 | country_code == 18 | 
                country_code == 19 | country_code == 21 | country_code == 22 | country_code == 24 | country_code == 26 | 
                country_code == 28 | country_code == 29 | country_code == 30 | country_code == 32 | country_code == 33 | 
                country_code == 34 | country_code == 35)

LAI_meanEU <- LAI_EU %>%
  summarise(mean = mean(LAI_Index_D7w_2020s))
LAI_meanEU$group <- "EU mean"

LAI_sample <- LAI_EU %>% filter(country_code == 7 | country_code == 8 | country_code == 11 | country_code == 13 | country_code == 14 | 
                           country_code == 15 | country_code == 18 | country_code == 28 | country_code == 34)

LAI_meansample <- LAI_sample %>%
  summarise(mean = mean(LAI_Index_D7w_2020s))
LAI_meansample$group <- "Sample mean"

LAI_meanfin <- rbind(LAI_meanEU, LAI_meansample)

FigA5 <- ggplot(data=LAI_EU, aes(x=reorder(Country_name_57, LAI_Index_D7w_2020s), y=LAI_Index_D7w_2020s)) +
         geom_bar(stat="identity")

FigA5 <- FigA5 + geom_hline(data=LAI_meanfin, aes(yintercept=mean, linetype=group))

FigA5 <- FigA5 + coord_flip() 

FigA5 <- FigA5 + ylab("Local Autonomy Index (LAI) in 2020") + xlab("") 

FigA5 <- FigA5 + scale_color_grey() + scale_fill_grey() + theme_own2

FigA5

# 7.4 Table A6 & Figure A6 ################################################
TableA6 <- list(
  "(1)" = lm(trstdiff ~ resent_munis.std*peri + cntry, data = dat_full),
  "(2)" = lm(trstdiff ~ resent_munis.std*peri + gender + age.std + educ + inc_dec + cntry, data = dat_full),
  "(3)" = lm(trstdiff ~ resent_munis.std*peri + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = dat_full),
  "(4)" = lm(trstdiff ~ resent_munis.std*peri + resent_munis.std*rur + cntry, data = dat_full),
  "(5)" = lm(trstdiff ~ resent_munis.std*peri + resent_munis.std*rur + gender + age.std + educ + inc_dec + cntry, data = dat_full),
  "(6)" = lm(trstdiff ~ resent_munis.std*peri + resent_munis.std*rur + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = dat_full),
  "(7)" = lm(trstdiff ~ resent_munis.std*peri*rur + cntry, data = dat_full),
  "(8)" = lm(trstdiff ~ resent_munis.std*peri*rur + gender + age.std + educ + inc_dec + cntry, data = dat_full),
  "(9)" = lm(trstdiff ~ resent_munis.std*peri*rur + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = dat_full)
)

FE_add3 = data.frame(matrix(c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), ncol = 10)) 
attr(FE_add3, "position") = 33

modelsummary(TableA6, fmt = 3, stars = T,
             title = 'OLS regression results: Local trust differential on place-based resentment scale conditional on urban-rural residence and centre-periphery residence.', 
             add_rows = FE_add3,
             coef_omit = "cntry", 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "periperiphery" = "Centre-periphery (b.=centre)",
                          "rurRural residence" = "Rural residence (b.=urban residence)",
                          "genderfemale" = "Gender (b.=male)",
                          "age.std" = "Age (Std.)",
                          "educHigh" = "Education (b.=low)",
                          "inc_dec" = "Income (Deciles)",
                          "lr.std" = "Left-right (Std.)",
                          "soctrst.std" = "Social trust (Std.)",
                          "plid.std" = "Place-based identity (Std.)",
                          "vote_oppopposition" = "Opposition vote (b.= government vote)",
                          "resent_munis.std:periperiphery" = "Place-based resentment (Std.) X Periphery residence",
                          "resent_munis.std:rurRural residence" = "Place-based resentment (Std.) X Rural residence", 
                          "periperiphery:rurRural residence" = "Periphery residence x Rural residence",
                          "resent_munis.std:periperiphery:rurRural residence" = "Place-based resentment (Std.) X Periphery residence x Rural residence",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

threewayint <- lm(trstdiff ~ resent_munis.std*peri*rur + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = dat_full)

FigA6 <- interact_plot(model = threewayint, pred = resent_munis.std, modx = rur, mod2 = peri, interval = TRUE, int.type = "confidence",
                      int.width = 0.95, colors = "Greys", mod2.labels = c("Centre", "Periphery")) 

FigA6 <- FigA6 + 
  labs(y = "Local trust differential",
       x = "Place-based resentment (Std.)") + theme_own2

FigA6 <- FigA6 + labs(x = "Place-based resentment (Std.)", y = "Local trust differential", title = "")

FigA6

# 7.5 Figure A7 ###########################################################
set.seed(270923)
desc_full <- dat_full

desc_rur <- desc_full %>% filter(rur == "Rural residence")
desc_rur$pbrter <- cut(rank(desc_rur$resent_munis.std, ties = "random"), 3, lab = FALSE)

desc_low <- desc_rur %>% filter(pbrter==1)
dat_trst <- desc_low %>% select(pers_id, nattrst, loctrst)

summary(dat_trst$nattrst)
summary(dat_trst$loctrst)

t.test(dat_trst$loctrst, dat_trst$nattrst, paired = TRUE, alternative = "two.sided") 

dat_trst_long <- melt(dat_trst, id=c("pers_id")) 
dat_trst_long$variable <- ifelse(dat_trst_long$variable=="nattrst", "National trust", "Local trust")

FigA7A <- ggplot(dat_trst_long, aes(variable, value)) +
  stat_boxplot(aes(variable, value), width=0.5)+  
  geom_boxplot(aes(variable, value),outlier.shape=1, width=0.5) +    
  stat_summary(fun=mean, geom="point", size=2, color="#FF0000") + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", color = "#FF0000", width = 0.1) + 
  scale_color_manual(values = c("#666666", "#000000")) +
  stat_compare_means(label = "p.format", 
                     method = "t.test", paired = TRUE, comparisons = list(c("National trust", "Local trust"))) + 
  labs(y = "Local trust / National trust", title = "A: Low resentment \n(1st tertile)") + theme_own1

desc_high <- desc_rur %>% filter(pbrter==3)
dat_trst <- desc_high %>% select(pers_id, nattrst, loctrst)

summary(dat_trst$nattrst)
summary(dat_trst$loctrst)

t.test(dat_trst$loctrst, dat_trst$nattrst, paired = TRUE, alternative = "two.sided") 

dat_trst_long <- melt(dat_trst, id=c("pers_id")) 
dat_trst_long$variable <- ifelse(dat_trst_long$variable=="nattrst", "National trust", "Local trust")

FigA7B <- ggplot(dat_trst_long, aes(variable, value)) +
  stat_boxplot(aes(variable, value), width=0.5)+  
  geom_boxplot(aes(variable, value),outlier.shape=1, width=0.5) +    
  stat_summary(fun=mean, geom="point", size=2, color="#FF0000") + 
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", color = "#FF0000", width = 0.1) + 
  scale_color_manual(values = c("#666666", "#000000")) +
  stat_compare_means(label = "p.format", 
                     method = "t.test", paired = TRUE, comparisons = list(c("National trust", "Local trust"))) + 
  labs(y = "Local trust / National trust", title = "B: High resentment \n(3rd tertile)") + theme_own1

loc_rur <- lm(loctrst ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
dat_loc_rur <- broom::tidy(loc_rur)
dat_loc_rur %<>% filter(term=="resent_munis.std")
dat_loc_rur$typeDV <- "Local trust"

nat_rur <- lm(nattrst ~ resent_munis.std + gender + age.std + educ + inc_dec + lr.std + soctrst.std + plid.std + vote_opp + cntry, data = subset(dat_full, dat_full$rur=="Rural residence"))
dat_nat_rur <- broom::tidy(nat_rur)
dat_nat_rur %<>% filter(term=="resent_munis.std")
dat_nat_rur$typeDV <- "National trust"

dat_comp_fin <- rbind(dat_loc_rur, dat_nat_rur)
dat_comp_fin$var <- "Place-based resentment (std.)"

dat_comp_fin$typeDV <- as.factor(dat_comp_fin$typeDV)
dat_comp_fin$estimate <- as.numeric(dat_comp_fin$estimate)
dat_comp_fin$std.error <- as.numeric(dat_comp_fin$std.error)

FigA7C <- ggplot(dat_comp_fin, aes(x = typeDV, y = estimate))
FigA7C <- FigA7C + geom_hline(yintercept = 0, color = "red", lty = 2)

FigA7C <- FigA7C + geom_linerange(aes(x = typeDV, ymin = estimate - std.error*interval1,
                                ymax = estimate + std.error*interval1),
                                linewidth = 1, position = position_dodge(width = -1/2))

FigA7C <- FigA7C + geom_pointrange(aes(x = typeDV, y = estimate, ymin = estimate - std.error*interval2,
                                 ymax = estimate + std.error*interval2, shape = typeDV),
                                 linewidth = 1/2, position = position_dodge(width = -1/2)) 

FigA7C <- FigA7C + xlab("DV: Local trust/National trust") + ylab("Estimated coefficient for place-based resentment (Std.)") + ggtitle("C") + coord_flip() + theme_bw(base_size = 20)

FigA7C <- FigA7C +  
  theme(legend.position="bottom") +
  theme(legend.title=element_blank())

ggarrange(ggarrange(FigA7A, FigA7B, ncol = 2),
          FigA7C,                                              
          nrow = 2)

# 7.6 Table A9 & Table A10 #################################################
TableA9 <- list(
  "(1)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="DE" & dat_full$rur=="Rural residence")),
  "(2)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="DE" & dat_full$rur=="Urban residence")),
  "(3)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="ES" & dat_full$rur=="Rural residence")),
  "(4)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="ES" & dat_full$rur=="Urban residence")),
  "(5)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="DK" & dat_full$rur=="Rural residence")),
  "(6)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="DK" & dat_full$rur=="Urban residence")),
  "(7)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="FR" & dat_full$rur=="Rural residence")),
  "(8)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="FR" & dat_full$rur=="Urban residence"))
)

modelsummary(TableA9, fmt = 3, stars = T,
             title = 'OLS regression results: Local trust differential on place-based resentment scale (urban vs rural; high local autonomy; country-level).', 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))

TableA10 <- list(
  "(1)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="HU" & dat_full$rur=="Rural residence")),
  "(2)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="HU" & dat_full$rur=="Urban residence")),
  "(3)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="CZ" & dat_full$rur=="Rural residence")),
  "(4)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="CZ" & dat_full$rur=="Urban residence")),
  "(5)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="PL" & dat_full$rur=="Rural residence")),
  "(6)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="PL" & dat_full$rur=="Urban residence")),
  "(7)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="EL" & dat_full$rur=="Rural residence")),
  "(8)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="EL" & dat_full$rur=="Urban residence")),
  "(9)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="IT" & dat_full$rur=="Rural residence")),
  "(10)" = lm(trstdiff ~ resent_munis.std, subset(dat_full, dat_full$cntry=="IT" & dat_full$rur=="Urban residence"))
)

modelsummary(TableA10, fmt = 3, stars = T,
             title = 'OLS regression results: Local trust differential on place-based resentment scale (urban vs rural; low local autonomy; country-level).', 
             coef_map = c("resent_munis.std" = "Place-based resentment (Std.)",
                          "(Intercept)" = "Constant"),
             gof_map = c("nobs", "r.squared", "adj.r.squared"))
